pysolve tutorial

Introduction

pysolve is a Python library written to aid in solving systems of (nonlinear) equations.

This tutorial will walk through and explain in greater detail how to use pysolve.

Installation

To install

pip install pysolve

Note that there are external dependencies. pysolve is dependent on the sympy and numpy packages.

Specifying the model

The sample simulates the SIM model from "Monetary Economics" by Godley and Lavoie (2006).

There are four main parts when creating a model.

• Creation of the model object
• Specifying the variables
• Specifying the parameters
• Specifying the equations

Model

Sample line 4 creates the model object. Currently, there are no parameters to the model object.

model = Model()

Model attributes

• variables
• parameters
• equations
• solutions

Variables

A variable, within pysolve, is a value that the solver is allowed to change when solving a system of equations. It also means that an equation defining the variable is required before attempting to solve the system.

The solver will attempt to "solve" the equation by iterating, through whatever method is selected, until the value of all variables are within a specified tolerance.

Variable specification is done in Sample lines 6-17.

model.set_var_default(0)
    model.var('Cd', desc='Consumption goods demand by households')
    model.var('Cs', desc='Consumption goods supply')
    model.var('Gs', desc='Government goods, supply')
    model.var('Hh', desc='Cash money held by households')
    model.var('Hs', desc='Cash money supplied by the government')
    model.var('Nd', desc='Demand for labor')
    model.var('Ns', desc='Supply of labor')
    model.var('Td', desc='Taxes, demand')
    model.var('Ts', desc='Taxes, supply')
    model.var('Y', desc='Income = GDP')
    model.var('YD', desc='Disposable income of households')

The first parameter (and only required parameter) is the name of the variable. This is case-sensitive and must be unique (across all variables and parameters) within the model. In general, avoid using common mathematical names, such as I, E, pi, and nan, as these are used by sympy.

Accessing the variables

The variable is returned from the call to model.var() and this can be used to modify the variable. In addition, the variables can be retrieved from the model also. The following snippets of code will return the same variable for x:

x_var = model.var('x')

is equivalent to

model.var('x')
    x_var = model.variables['x']

variable Attributes

• name
• desc
• default
• value
• model
• equation

Defaults

It is highly recommended, though not required, that defaults be given to the variables. This is done on Sample line 6.

model.set_var_default(0)

The default value is used to give the variable an initial value when the variable is created. If variables are created before the call to _set_vardefault(), they will not have an initial value. If the variable is used before a value is assigned, this may cause hard-to-debug errors when running the simulation.

Default values are also used in cases where values for a previous iteration are requested (but we may not have iterated far enough).

Defaults may also be given to individual variables using the default named parameter. This value will take precedence over the _set_vardefault value.

model.var('YD', desc='Disposable income', default=1.2)

Creating multiple variables at once

The creation of multiple variables at the same is also supported. However the ability to specify descriptions and defaults is not allowed. Although those values can be set separately.

model.vars('x', 'y', 'z')

    # Note that setting the default will not set the value for x, as that is
    # done only during creation of the variable.
    model.variables['x'].default = 1.1
    model.variables['x'].description = 'variable for x'
    model.variables['x'].value = 2

Parameters

A parameter, within pysolve, is a value that the solver is not allowed to change when solving a system of equations. That is, it is considered a constant by the solver. An error will be generated if an equation is used to define the value of a parameter. Note that I use the term "parameter" for true system parameters and exogenous variables.

Note that the solver will not change the value of a parameter. But the value of a parameter may be changed in between iterations by the user.

Parameter specification is done in Sample lines 19-23.

model.param('Gd', desc='Government goods, demand')
    model.param('W', desc='Wage rate')
    model.param('alpha1', desc='Propensity to consume out of income')
    model.param('alpha2', desc='Propensity to consume out of wealth')
    model.param('theta', desc='Tax rate')

The first parameter (and only required parameter) is the name of the parameter. This is case-sensitive and must be unique (across all variables and parameters) within the model. In general, avoid using common mathematical names, such as I, E, pi, and nan, as these are used by sympy.

Accessing the parameters

The parameter is returned from the call to model.param() and this can be used to modify the parameter. In addition, the parameters can be retrieved from the model also. The following snippets of code will return the same parameter for a:

a_param = model.var('a')

is equivalent to

a_param = model.variables['a']

parameter Attributes

• name
• desc
• default
• value
• model

Defaults

It is highly recommended, though not required, that defaults be given to the parameters. For example,

model.set_param_default(0)

The default value is used to give the parameter an initial value when the variable is created. If parameters are created before the call to _set_paramdefault(), they will not have an initial value. If the parameter is used before a value is assigned, this may cause hard-to-debug errors when running the simulation.

Default values are also used in cases where values for a previous iteration are requested (but we may not have iterated far enough).

Defaults may also be given to individual parameters using the default named parameter. This value will take precedence over the _set_paramdefault value.

model.param('YD', desc='Disposable income', default=1.2)

Equations

The equations are at the heart of the model. However there are some specifics that must be understood when writing the equations. The equations may be non-linear, there is no linearity requirement. Because of this, the equations may or may not converge to a solution.

Hers are the equations from Sample lines 25-35

model.add('Cs = Cd')  # 3.1
    model.add('Gs = Gd')  # 3.2
    model.add('Ts = Td')  # 3.3
    model.add('Ns = Nd')  # 3.4
    model.add('YD = (W*Ns) - Ts') # 3.5
    model.add('Td = theta * W * Ns')  # 3.6, theta < 1.0
    model.add('Cd = alpha1*YD + alpha2*Hh(-1)') # 3.7, 0 < alpha2 < alpha1 < 1
    model.add('Hs - Hs(-1) =  Gd - Td')  # 3.8
    model.add('Hh - Hh(-1) = YD - Cd') # 3.9
    model.add('Y = Cs + Gs') # 3.10
    model.add('Nd = Y/W') # 3.11

The model equations follow the python way of writing equations. Thus '*' is used for multiplication for example.

Requirements

• Each variable can only be solved for once

  There cannot be two equations solving for a variable. For example, the following is not allowed

  x = y + z + ...
    x = log(y) + 2*z + ...
  

• Each variable must have an equation defining it

  There must be an equation that "defines" the variable. This is mostly a requirement for the Gauss-Seidel algorithm. Thus, for a variable x, there must be an equation of the form, x = ...

• The left-hand side of the equation must contain only one variable

  The variable on the left-hand side of the equation is the variable being defined. Thus there can only be one variable on the left-hand side.

  This also means that there can be constants and parameters on the left-hand side. For example:

# Good
    10*x = .....
    x - x(-1) = ....

    # Bad
    x*y = ....
    log(x) = ...

• The variable on the left-hand side must be of linear form

  There cannot be an exponential power or inside of a function parameter.

# Bad
    log(x) = ...
    x**2 = ...

• Only variables can be solved for (no parameters)

  Only variables can be solved for. Parameters cannot be defined.

Accessing previous iterations's values

In the sample, there are some lines which look like this:

model.add('Hh - Hh(-1) = YD - Cd') # 3.9

This equation uses the value of Hh from a previous iteration (or solution from a previous iteration). This is done with the "Hh(-1)". Using "-1" will refer to the previous iteration, "-2" will refer to the iteration before that, and so on.

Values from a previous iteration are treated as parameters by the system.

This will work for both parameters and variables.

Functions available for use within an equation

• The delta ('d') function

  This is a special function. The purpose of this is to reduce the amount of work in order to display the difference between the current value and the immediate past value.

  Thus

  d(x)
  

  is equivalent to using

  x - x(-1)
  

  However, "d(x)" cannot be used on the left-hand side of an equation.

• if_true

  This function will return the value 1.0 if the argument is true, 0.0 otherwise.

  x = 5
    if_true(x >= 4)
  

  The call to if_true() here will return 1.0.

• Abs

  Absolute value

• exp

  The natural exponential function.

• log

  The natural logarithim.

• Max

  Returns the maximum of two values.

• Min

  Returns the mininum of two values.

• sqrt

  Returns the square root of a value.

• sign

  Returns the sign of the argument.

Model.evaluate()

This is a special function that is available. This takes a string argument. It will evaluate the string and then return the value.

# assuming that x = 1 and y = 2
    # this will return the value 11
    val = model.evaluate('3*x + 4*y')

Equation attributes

• equation
• desc
• model
• expr
• func
• variable

Solving the model

To solve the model, we will run various iterative algorithms until the values converge. In addition, the Model object will store solutions from previous runs.

Sample

Setting initial values

Use the _setvalues() function to change the values for a group of parameters for variables. Arguments to this function can either be a python dictionary() or a python list(). There are two use cases here.

The specification of the value involves the name of the variable and the value. The value may either be a numeric value (integer or float) or may be a string value. In the case of a string value, the string will be evaluated with the current context of values.

The python dictionary is the more natural, however, python dictionaries are unordered, so the actual order in which the variables are applied is not apparent.

If a list ot tuples is passed in, the values will be evaluated in the order of the list, thus one can use the value of previous entries and use calculated values.

Example using a dictionary

model.set_values({'x': 1.1,
                      'y': 2.2})

Example using a list of tuples

model.set_values([('x', 1.1),
                      ('y', 'x*2')]

Using a formula for y in the dictionary example is not guaranteed to work since the order of the dictionary is not guaranteed.

Invoking solve()

model.solve(iterations=10, threshold=0.001, method='gauss-seidel', debuglist=None)

The method is pretty straightforward.

• iterations

  This is the maximum number of iterations allowed until we reach convergence. If all of the variables fail to converge, then a SolutionNotFound() exception is raised.

• threshold

  This is the relative tolerance that is used to test for convergence.

• method

  This is the method used to find a solution. There are currently three options: 'gauss-seidel', 'newton-raphson', and 'broyden'. For 'gauss-seidel', the equations are evaluated in the order in which they are added.

• debuglist

  If a list is passed in, then the debuglist will contain a list of the values that the solver is using, one entry per iteration.

Using the model

Using the model refers to being able to access the solutions.

Sample

Here is some code that prints out the values from the zeroth iteration (the initial values), the first to third iterations, and the final iteration.

Note that the values are not quite whole. So at the end, I take the final solution, round off the values to 1 decimal place and then print it out. Note that _printsolutions() expects a list, which is why I convert _nicesolution into a list before invoking the function.

The solutions data

An individual solution is just a python dictionary, mapping the variable/parameter name to a value.

Accessing previous solutions

The data in model.solutions is a list of solutions, where each individual solution is a dict() of parameter and variable values.

Thus, the normal python list access and slicing features are available to access the information.

So to access the last solution in the list

model.solutions[-1]

To gather solutions from 5 to 10 (python lists are 0-based)

model.solutions[5:10]

To get all solutions starting from 5 on

model.solutions[5:]